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1.0 Introduction 



The problem which we will be considering is that of interpolation in 
two (or more) variables, where the data to be interpolated is at irregularly 
spaced points. While the discussion will often center on functions of two 
variables, most of the ideas are readily extendable to any number of inde- 



pendent variables. 

Assume that data points (X^, f^) , k = 1,...,N are specified, where 
the lie in n-space. We make no assumption about the disposition of 

the points X^ . We wish to construct a function F(X) which is defined 
in some region which contains all the X^ and such that F(X^) “ f^ , 
k = 1,...,N . We wish to discuss a broad enough class of interpolating 
functions so that we can require that F(X) is not only continuous, but 
perhaps has some continuous derivatives as well. 

A number of methods have been proposed for the problem. For the most 
part these are global methods, most of which require the solution of a sys- 
tem of N linear equations for the coefficients of a set of basis functions. 
Typical of these are polynomial interpolation [12], and optimal approximations 
in various spaces of functions [6], [7], [9], [2]. 

Another type of global approximation (although the computational version 

is local) is the idea of taking the function value to be a weighted average 

of the values at data points, the weight being a decreasing function of dis- 

-2 

tance [11]. A typical weight function is w(d) = d . There are several 
serious faults with this approximation, such as the first partial derivative 
being zero at each data point. These faults are overcome by the imposition 
of somewhat artificial conditions on a variation of the approximation. A 
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remaining bad feature is that the approximation is dependent on the coordi- 
nate system, or rather, the measure of distance. 

A variation of the distance weighting idea is proposed by McLain [5]. 
There the idea is to fit the data with an equation of specified type in the 
weighted least squares sense. The weight attached to each point is a de- 
creasing function of its distance from the point at which the approximation 
is desired. McLain investigated a number of different fitting functions and 
weight functions. This approach eliminates in a natural way most of the bad 
features of Shepard’s approximation, although not the dependence on the co- 
ordinate system. McLain claims the approximation has partial derivatives of 
all orders, and although this seems plausible, it is neither obvious nor eas- 
ily verifiable. 

Where the data is on a somewhat regular set of points, e.g., if the grid 
is topologically equivalent to a regular grid [3], or if the grid can be tri- 
angulated, or represented as a union of convex quadrilaterals [4], special 
approaches may give good results. Approximations in these cases would gen- 
erally be local with limited smoothness, although Ferguson’s [3] parametric 
representation is global and a type of spline. 

Another approach to the interpolation problem, the idea of which is the 
basis for most of our investigation, is due to Maude [8]. The resulting ap- 
proximation can be made to have arbitrarily high smoothness and is locally 
determined. Since we will discuss the idea at length in the next section, 
it is mentioned here only for completeness. 

The ideas which will be important to our discussion will include the 
following: (i) No assumption will be made about the disposition of the data 

points, and in particular we do not preclude the possibility of a regular 
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grid of points; (ii) A smooth approximation in the sense that the function 
and some of its partial derivatives are continuous, is desirable; (iii) The 
number of points, N , may be so large that the approximation must be lo- 
cally determined, or at least have a local basis (such as the B-splines are 
a local basis for univariate spline approximation) to be computationally 
viable; (iv) It is desirable for the approximation to be independent of 
translation and contraction or expansion of the coordinate system; and (v) 

It is desirable that the procedure be robust so that failure of the exis- 
tence of the approximation can perhaps be circumvented by some slight al- 
teration. 

It is likely that (iii) rules out the optimal approximations because 
they are global, and the basis functions are sufficiently complex that it 
is quite unlikely a set of basis functions can be derived which have com- 
pact support. A large number of points also rules out polynomial approxi- 
mation. 

2.0 A Class of Interpolation Methods 

We will begin this section with a brief review of a method due to Maude 
[8]. Some of the difficulties experienced with the method will be mentioned. 
A general class of methods which contains Maude’s will be developed and some 
properties of the class of methods will be discussed. 

2. 1 Maude’s Method 

) and f = f 
k 

each point (x , y ) associate a circle C , of radius r with center at 
n’ n n 

(x , y ) . This circle is chosen so that the closed disk contains (x ,y ) 
n n n n 

and its five nearest neighbors from the set {(x^, y^) } , with no more than 
four of the neighbors in the interior. There may be ’ties’ for the position 



(x^, y^) , k = 1,...,N be specified. With 
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of fifth nearest point, and one needs some means of deciding which point on 
the boundary is to be used in the subsequent calculations. Let the second 
degree polynomial interpolating the function f(x,y) at (x^y^) and ^ ts 
five nearest neighbors be denoted by Q n (x,y) . Let 

| l-3s 2 + 2s 3 , 0 < s < 1 

w(s) ■ { 

| 0 , s > 1 



To obtain the value of the interpolating function F(x,y) at the point 

(x ,y ), we calculate a weighted average of the values of Q^(x ) for 

* * 

those values of n such that (x ,y ) lies in the circle . In particu- 

lar, 

f(x*,y*) = l w(V(x*-x n ) 2 + (y*-y n ) 2 /r n )Q n ^ x *» y *)/ 1 W (^ x * _x n ) 2 + (y*-y n ) 2 /r n 

n=l / n=l ' 



) 



Because 



* * 



is zero unless (x ,y ) lies inside 



the circle , the sum can be taken over only those values of n such 

* * * * 
that (x ,y ) is in . The function is defined for all points (x ,y ) 

which lie in the interior of one or more of the circles. By virtue of the 

weights being zero, with zero first partial derivatives, on the boundary of 

the circles, the interpolating function F is continuous, with continuous 

first partial derivatives. 

We can note that F is a bivariate analogue of the quintic spline of 
deficiency three [1]. Unfortunately, unlike the univariate case, the bivar- 
iate version above can fail. It is not necessarily true that the polynomials 
must exist. In fact, on a uniform grid a boundary point and its five 
nearest neighbors will lie on two straight lines, hence the interpolating 
polynomial does not generally exist. Non-existence of a particular can 

be bypassed by simply ignoring that circle in the definition of F . This is 
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handled computationally by setting = 0 • For isolated instances the effect 

may be minimal, although failure for very many points could severely curtail 
the region of definition of the interpolating function and even exclude some 
data points. 

Experience with this version of the computational scheme indicates that 
the behavior of the are adversely affected by the presence of even 

moderately steep gradients. This problem is not unique to polynomials, but 
they do seem to be more sensitive than other basis functions we have investi- 
gated. In addition, regions where the data points are relatively sparse seem 
to yield poor approximations. 

2. 2 A Generalization 

We will develop a class of interpolation methods based on Maude 1 s idea, 
which will contain his method as a special case. We will begin with a weight- 
ing function associated with each point, although in the application it will 
probably not be the starting point. 

We presume that data points corresponding function values 

f, , k = 1 N have been given. With each point (x ,y ) associate a 

k n n 

nonnegative weight function w^(x,y) which has compact support S n , such 

that (x ,y ) is contained in the interior of S • Let U (x,y) be a local 
n n n n 

interpolating function which is defined on and interpolates to the value 

f fc at each ( x k »y fc ) in . 

We define the value of the interpolating function F(x,y) at the point 
& & 

(x ,y ) to be 

N N 

* * r * * * * „ * * 

f(x ,y ) = 2, w n ( x »y )u n (x »y ) / I w n ( x >y ) » 

n=l n=l 
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provided (x ,y ) lies in the interior of one of the . Since pre- 

sumably (x ,y ) does not lie in all , we may want to write 



, * * 
F(x ,y ) 



l w (x ,y )U (x ,y )/ 2 , W (x ,y ) 
nel* nel* 



where I = { n : (x ,y ) e interior (S ) } . 

We can make a number of observations about the function F . 

(a) FCxj^.y^ = f fc since if (x^y^ € interior of S n , u n (\»y k ) = f k 5 

(b) F is at least as smooth as the least smooth of the functions w^ and 

* * 



, k = 1,...,N ; (c) F(x ,y ) is locally determined by the points 

(x_ ,y. ) e u S ; (d) the error in the interpolation is no worse than 
k k T i n 

nel* 

that of the poorest interpolation by a , since the value of F is a 

convex combination of the , hence the error is a convex combination of 

the local interpolating errors* 

It is trivially observed that Maude’s method is a special case, pro- 



vided the (U^ above) exist. For Maude’s method the error is of the 

2 

form 0(h ) , where h = max r. [12]* If one considers a sequence of ap- 

k k 



proximations as more points are added, in a fashion such that h -> 0 , it 

is necessary that the ’condition’ of the sets of nearest neighbors not de- 

2 

teriorate too much in order to maintain the 0(h) error estimate* 



3.0 Implementation 

In order for a scheme such as that outlined in the previous section to 
be reasonably easy to implement, the weight functions, in connection with 
the shape of the support regions must be chosen with some care. Since 

we wish to be able to incorporate smoothness, which requires the value of 
w^ and some of its partial derivatives to be zero on the boundary of , 

we are led rather naturally to such uncomplicated shapes as squares and 
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circles. Squares are associated with and norms on 2-space and 

circles with the norm. The weight functions for these regions can be 

easily built up from univariate functions, as follows. Let G3(s) denote a 
univariate function with appropriate smoothness. For example 



0 )( S ) 



l-|s| 

0 



|s| 1 1 

I s I > 1 



is continuous, while 



«(s) 



l-3s^ + 2 1 s | | s | <_ 1 

0 | s | > 1 



has continuous derivatives as well. Similar functions with arbitrarily 
high smoothness can be generated. Reasonable weight functions for balls 
of radius r^ centered at (x^y^ are given by! 







/ (X-X k } + 


(y-y k ) \ 


h • 


^(x.y) = 


0)1 

V r k 




i 


w k (x,y) = 


“ (y ( - x ~\ )2 


+ 


l : 
00 


w k < x > y ) = 


3 

1 u 
X 

3 


m 



- (y-y k ) 



The choice of region will affect the approximation slightly. Which should 
be used is a matter of personal preference, in the author's opinion. The 
square regions have the distinction that the weight functions are piecewise 
polynomials, which results in the function F being a rational function if 
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the are polynomials or rational functions. While this is aestheti- 

cally more pleasing than a function involving square roots as is obtained 
for circles, the author believes the practical difference to be minimal. 

It may be an advantage to use regions which do not have flat sides. Some 
calculations have been done both ways with neither one seeming to be sig- 
nificantly better than the other. 

We will generally assume that the regions are chosen to contain 

fc ) . This 

will allow the use of the same type of local interpolating function for 
each region and will simplify the coding of programs to implement the pro- 
cedure. The value of m will greatly influence the computational effort 
required, since the larger m is, the more the support regions S, will 

K. 

overlap. On the other hand, too small a value for m will not allow the 
local interpolating function to mimic f(x,y) suitably, and may cause holes 
in the region on which F(x,y) is defined, particularly if the data is 
sparsely located in some areas. Hence a balance must be struck between 
computational efficiency and adequate sampling of f(x,y) along with cov- 
erage considerations. 

The choice of local interpolating functions is very broad, since 

essentially the only restrictions are appropriate smoothness and interpola- 
tion of the data. The obvious choices are polynomials consisting of m 
monomials, chosen in some fashion. While this works well for some sets of 
data, it has not yielded the best approximations overall. 

Better results have been obtained by using the optimal approximations 
previously mentioned [6], [9], [2]. While optimal approximations have not 
received great attention in applications, they do have very desirable proper- 
ties and prove to yield quite good results, with one exception, when applied 



a set number of nearest neighbors, say m-1 , of the point (x^y 
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to our test problems in the next section. These functions are splines and 
thus are better able to adapt to certain data. 

When using optimal local interpolating functions it is probably de- 
sirable to use weight functions of no more smoothness than the interpolating 
functions. Thus if the interpolating function is the optimal approximation 
from the class ^ (see [10] for the definition of this space) one should 



choose 0)(s) 



1- I s I 

0 



I s I <_ 1 

I s I > 1 



since functions in 



B 



n,n 



are only continuous. 



4.0 Numerical Experiments 

It was deemed desirable to test a number of methods on a variety of sur- 
faces for purposes of comparison. Five sets of data were fit using eleven 
fitting functions, and one set of data was fit using ten of these procedures, 
the eleventh not being suitable. We will describe the procedures and the data 
sets and give a table of results with some discussion and observations. 

4. 1 The Data Sets 

Because of the variety of surfaces used by McLain [5], it was decided to 
use them for our tests also. Unlike McLain, we did not specify uniformly spaced 
points on [1,10] x [1,10] . Instead, in each of the 81 squares [i,i+l]x [j , j+1] 
i,j = 1,...,9 , a point was generated by a random number generator. These 

points are shown graphically in Figure 1. In addition, nine points were speci- 
fied at positions which were judged a priori to be somewhat critical points in 
defining the surface. Thus a total of 90 points were specified on surfaces 
S1-S5 on the square [1,10] x [1,10]. 
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Figure 1: Interpolation Points 
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We point out that McLain was interested in contours for topographical 
data and he describes the surfaces in topographical terms. We shall do this 
also as this leads to an easy visualization of the surfaces. 

(SL) This surface is a section of a sphere and is described by 

V 2 2 

64 - (x - 5.5) - (y - 5.5) . The additional points speci- 

fied were (a, 8) , a, 8 = 1* 5.5, 10. 

(S2) This surface is a hill rising steeply from a plain, and is described 

2 2 

by f(x,y) = exp(-(x-5) - (y-5) ) . The additional points specified were 

(4,5) , (4.5, 4.5) , (4.5, 5.5) , (5,4) , (5,5) , (5,6) , (5.5, 4.5) , (5.5, 5.5), 
(6,5) . 

(S3) This surface is similar to S2, a hill rising less steeply from a 

plain. It is described by f(x,y) = exp ^ j # The additional 

points specified were the same as for S2. 

(S4) This surface is a long narrow hill rising from a plain, and is de- 

2 2 

scribed by f(x,y) = exp(-(x+y - 11) - (x-y) /10) . The additional points 

were (3,8), (4.5, 6.5), (4.5, 5.5), (5.5, 4.5), (5.5, 5.5), (5.5, 6.5), 

(6.5, 4.5), (6.5, 5.5), (8,3) . 

(S5) This surface is a plain and a plateau separated by a sharp rise, 
almost a cliff. It is described by f(x,y) = tanh(x + y - 11) . The addi- 
tional points were (1,10), (2.5, 8.5), (4,7), (5,5), (5.5, 5.5), (6,6), (7,4), 
(8.5, 2.5), (10,1) . 

(S6) This surface is that given by Ferguson [3] and is not described 
mathematically. For purposes of better point spacing the y coordinate was 
multiplied by three. This was almost necessary because a point and its 
nearest five neighbors would often include five points very nearly on a 
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straight line when the original data was used. Multiplication of the y-coordin* 
by a constant is equivalent to using ellipsoidal regions in methods Fl - F5. II 
has no effect on methods F6 - F8, and has a very beneficial effect on methods 
F9 - F10- The data is listed in Table 1. 

4. 2 The Interpolation Functions 

The eleven interpolation functions consist of five different implemen- 
tations of the method described in Section 2, three global methods for 
comparison purposes, and three versions of McLain’s distance weighted least 
squares approximation. We will discuss each of these in turn. Several of 
them involve use of optimal approximations in certain spaces of functions, and 
references to more information on these topics are [10], [6], [9], and [2]. 

(Fl) This is Maude’s method, using as a second degree polynomial 

interpolating to the point and its five nearest neighbors. This was described 
in Section 2.1. 

(F2) This is a method of Section 2, using the optimal approximation 
from the class of functions ^ for the local interpolating function for 

a point and its five nearest neighbors. Six points per region was chosen as 
being sufficient for coverage, and adequate for function definition, as well 
as comparable computationally to Maude’s method. Circular regions were used 



Since Barnhill and Nielson [2] have shown that Sard spaces of type B* 
have a reproducing kernel, one can find the optimal approximation as a linear 
combination of basis functions (or representers of the point evaluation 
functionals), one of which is associated with each data point. The basis 




Is! < i 



with w(s) 



0 



|s| > 1 
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X 


y 


f 


2.0 


15.0 


2.0 


2.49 


7.647 


2.7 


2.981 


0.291 


2.9 


3.471 


- 7.062 


3.0 


3.961 


- 14.418 


3.0 


7.45 


12.003 


2.0 


7.35 


6.012 


3.0 


7.251 


0.018 


2.5 


7.151 


- 5.973 


1.5 


7.051 


- 11.967 


2.0 


10.901 


9.015 


1.5 


10.751 


4.536 


1.425 


10.602 


.06 


1.35 


10.453 


- 4.419 


1.276 


10.304 


- 8.895 


1.2 


14.055 


10.509 


1.0 


14.194 


6.783 


0.8 


14.331 


3.054 


1.2 


14.469 


0.672 


1.6 


14.607 


- 4.398 


1.25 


15.0 


12.0 


0.0 


15.729 


8.067 


0.0 


16.457 


4.134 


0 .20 


17.185 


0.198 


0 .60 


17.914 


- 3.735 


1.2 



Table 1 
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function for the point (x^.y^) in n is 

^(x.y) = [l + (x-a) - ( x -x k ) + ][i + (y-b) - (y-y k > + 3 . 



where 




0 , t £ 0 

. is the truncated power function and 

t , t > 0 



(a,b) is a parameter pair which in this case must satisfy >_ a , 

y > b for n = 1,2,...,N . In general (a,b) need not satisfy this 
n 

condition, but we have simplified by requiring the condition. 

We can note that optimal approximations from B n,n are piecewise 
bilinear functions with knots along lines through the data points and paral- 
lel to the axes. 

(F3) This is a method of Section 2, similar to function F2, except that 
is taken to be the optimal approximation from 21 9 anc ^ we: *-8kt 

I I- 3s 2 + 2 1 s | | s | £ 1 

0 | s | > 1 



Again, the nearest five neighbors in a circular region were used. The basis 
function associated with the point (x^,y^) in this case is 



where 



K^x.y) = g(a,x k ,x)g(b,y k ,y) 

(x^-a) (x-a) 2 

g (a.x^.x) = 1 + (x^ - a) (x-a) + 



(x-a) 3 (x-Xj^) 3 

”31 + 31 



where the notation and parameters are as described before. 

(F4) This approximation is the same as F3 except that the nine nearest 
neighbors in a circular region were used to define the local interpolating 
function. This is to test the effect of including more points in the local 
approximation. 
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(F5) This is a method of Section 2, where the is taken to be the 

optimal approximation from the space of functions described by 

Mansfield [6]. Again, a circular region containing a point and its five 

I i- | S | , I s I £ 1 

nearest neighbors was used, with w(s) = { . This 

(0 , |s| > 1 



space is constrained so as to give an exact approximation if the approxi- 
mated function is linear. Our version used the basis functions given by 
Nielson [9] and are of the form 






( *'V+ 

3! 



+ 



(y - y k } + 

3! 



+ 



[ ( \ - x) +" (a - x) + - ( \ ■ a) +] [ ( v y) + ■ (b ' y) + • <y k' b> +] * 

where a and b are parameters. The optimal approximation is a linear 

combination of and the . 

(F6) This is the global optimal approximation from B u,n * 

(F7) This is the global optimal approximation from ® p 2 21 * 

(F8) This is the global optimal approximation from T^*^ . 

(F9) This is McLain's method M10 which he describes as being accurate 

with a satisfactory amount of computation. The weight attached to a data 

2 2 2 2 2 

point (x^.y^) w (d) = exp(d )/d where d = (x-x^) + (y-y^) is the 

distance squared from the point (x,y) at which the function value is desired. 
The fitting function is a second degree polynomial. 

(F10) This is the same as F9 except that distance is changed. Here 
2 2 2 

d = [ (x-x^) + (y-y k ) ]/10 . This is equivalent to shrinking the coordinates 
by a factor VlO , and is to test the effect of a different coordinate system. 
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(Fll) This is the same as F9 except that distance is changed. Here 



2 f 2 21 

d =10 (x-x^) + (y-y ) . This is equivalent to expanding the coordi- 
nates by a factor Ho and is to test the effect of a different coordinate 
system. 



4.3 Comparisons 

For surfaces S1-S5 and all eleven interpolating functions the devi- 
ation at the 361 points (x,y) , x,y = 1(.5)10 was generated. The maximum 
deviation, the mean deviation, and the root-mean-square deviations are listed 
in Table 2. The number in parenthesis below each function is the approximate 
calculation time in seconds for all calculated data points. The calculations 
were done on the IBM 360/67 at the Naval Postgraduate School. These times 
are given only as an indication of the number of computations required for 
the various methods. 

Generally we can observe several things from the table: 

(i) Maude’s method (FI) does not compare favorably with any of the 
other methods, except when applied to surface SI . 

(ii) Surface SI seems to be relatively difficult to fit, although 
Maude’s method (FI) and McLain’s method (F9) work very well. 

(iii) The local weighted optimal approximations (F2-F5) often are more 
accurate than the corresponding global versions (F6-F8) . A notable excep- 
tion is the application to surface SI . 

(iv) McLain’s method (F9) does very well on all these surfaces and 
is often the most accurate. The suspicion that the distance weighting func- 
tion is rather finely tuned for this particular coordinate system is rein- 
forced by poorer performance of the expanded and contracted length versions 
of it (F10-F11). The amount of computation required is considerably greater 
for the method than for any other methods considered here, even the global 
methods . 
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Surface 


SI 


S2 


S3 


S4 


S5 


Fitting Function 












FI : Maude 
(34) 


.06562 

.00390 

.00791 


1.390 

.06882 

.2106 


1.198 

.03964 

.1262 


3.798 

.09371 

.3789 


.7438 
. 04440 
.10468 


F2 : Bj.^ local 

(51) 


1.768 

.1378 

.2456 


.1704 

.00943 

.02661 


.1520 

.01180 

.02429 


.3393 

.01732 

.04202 


.4926 

.04828 

.09124 


F3 : B r2 21 local 
(63) 


1.284 

.1049 

.1859 


.1390 

.00726 

.01993 


.08185 

.00832 

.01483 


.2701 

.01446 

.03676 


1.111 

.03876 

.09576 


F4 = B f2,21 l0cal 
(143) 


.3472 

.03970 

.06454 


.1315 

.00727 

.01953 


.05723 

.00593 

.01073 


.2585 

.01758 

.03831 


1.178 

.04932 

.11336 


F5 : T 1 ' 1 local 
(60) 


.1068 

.01667 

.02477 


.1740 

.00884 

.02149 


.1143 

.01035 

.01797 


.2398 

.01493 

.03280 


.4656 

.03872 

.07300 


F6 : Bj.^ global 

(128) 


.2209 

.02724 

.05102 


.2829 

.05706 

.08218 


.1514 

.02234 

.03416 


.3547 

.05453 

.08126 


.5098 

.10322 

.14568 


F7 : B r2>2l global 
(154) 


.09152 

.00832 

.01465 


.4838 

.07447 

.1086 


.09318 

.01012 

.01719 


1.335 

.2005 

.3127 


.8164 

.10104 

.15686 


F8 : T 1,1 global 
(138) 


.1222 

.01355 

.02394 


.2693 

.05135 

.07601 


.1457 

.02001 

.03150 


.3549 

.05106 

.07797 


.4868 

.09632 

.13595 


F9 : McLain 
(372) 


.05223 

.00354 

.00827 


.1372 

.00824 

.02090 


.04894 

.00530 

.00873 


.1803 

.01124 

.02546 


.4164 

.03006 

.05960 


F10: McLain, d 2 «• d 2 /10 
(465) 


.04343 

.00621 

.00982 


.2158 

.02196 

.04363 


.08728 

.02174 

.03117 


.3894 

.02549 

.05174 


.2738 

.05374 

.07886 


2 2 

Fll: McLain, d ^Qd 

(87) 


5.735 

.03785 

.3169 


1.304 

.01222 

.07502 


.1254 

.00982 

.02010 


.4145 

.01539 

.03949 


4.118 

.04584 

.2318 



Table 2: Maximum, Mean and RMS deviations 
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(v) While functions in 
pi * 1 



71,11 



and T 



1.1 



are of similar smoothness 



(only continuous), T J ‘ >J ‘ gives considerably more flexibility which shows 
up when fitting surface SI. Results are comparable for surfaces S2-S5. 



(vi) Changing the number of neighbors included in a region seems to 
generally have only a small effect (F3 vs. F4). The one exception, again, 
is on surface SI . 

For surface S6 we could not compare exact errors since only data 
points were given. The data is apparently difficult to fit since it in- 
volves three relative maxima with a saddle point between two of them and 
a shallow valley between two of them, the entire surface only having 25 
points given on it. In this situation it is not surprising that interpo- 
lating surfaces may involve some overshoot. We list the minimum and maxi- 
mum observed values for functions F1-F10 in Table 3. Again, the flexi- 
bility of the space T^ > ^ is noted, and the fine tuning of the distance 
function for F9 , yielding a good result in FICh 
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Function 



Extreme Values 



FI 

F2 

F3 

F4 

F5 

F6 

F7 

F8 

F9 

F10 



-6.43, 7.18 
0 , 4.65 

0 , 5.52 

0 , 6.02 

0 , 3.00 

0 , 3.77 

0 , 5.05 

0 , 3.05 

-26.31, 25.83 
0 ,3.30 



Table 3 



5.0 Conclusion 



The main conclusion which can be made is that the method of Section 2 
yields good approximations provided the local interpolating functions are 
reasonably accurate. Six points per region appears to give sufficient func- 
tion definition and coverage at a reasonable computational load. Any of the 
optimal approximations seem to be satisfactory as local interpolating func- 
tions, and the use of them is strongly recommended over polynomials. 

The use of McLain’s method (F9) or a variation of it gives very good 
results if the computational load is not burdensome. Some tuning of the 
measure of distance may be necessary and is a definite disadvantage of the 
method. 

It appears that on certain surfaces, some methods will do poorly while 
on others the situation could be reversed. In terms of which of the methods 
discussed can ultimately be tuned to a given set of data, the author feels 
McLain’s method has a good capability. Unfortunately in case only isolated 
points are known the procedure for this tuning is not known, and it is proba- 
bly safer to use a method which adapts naturally to changes in coordinate 
systems and possible varying sparseness of data points, such as the methods 
of Section 2, where changing the obvious parameter (the number of neighbors 
included in a region) seems to have a fairly small influence. 
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